Systematic Computation of the Least Unstable Periodic Orbits in Chaotic Attractors 
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We show that a recently proposed numerical technique for the calculation of unstable periodic 
orbits in chaotic attractors is capable of finding the least unstable periodic orbits of any given order. 
This is achieved by introducing a modified dynamical system which has the same set of periodic 
orbits as the original chaotic system, but with a tuning parameter which is used to stabilize the orbits 
selectively. This technique is central for calculations using the stability criterion for the truncation of 
cycle expansions, which provide highly improved convergence of calculations of dynamical averages 
in generic chaotic attractors. The approach is demonstrated for the Henon attractor. 
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Unstable periodic orbits in chaotic attractors provide 
a useful hierarchical framework for calculations of prop- 
erties such as Lyapunov exponents, fractal dimensions 
and entropies of the attractors [Qj2) . Periodic orbits have 
been used to characterize the attractors in a variety of 
low dimensional dissipative dynamical systems, includ- 
ing model systems such as discrete maps as well as 
experimental time series In chaotic Hamiltonian 

systems, series expansions over unstable periodic orbits, 
within the semi-classical approximation, have been used 
to calculate the quantum energy level density as well as 
properties of the wave functions |t]]. Cycle expansion 
techniques gave rise to highly improved convergence, par- 
ticularly for systems in which the symbolic dynamics is 
well understood and long periodic orbits are well shad- 
owed by short ones 

The calculation of unstable periodic orbits in chaotic 
dynamical systems is a difficult computational problem. 
The difficulty is that in a chaotic system the numerical 
error grows exponentially with the length of the orbit. 
Therefore, only short unstable periodic orbits can be cal- 
culated using the standard techniques of map iteration. 
Moreover, even if some orbits of a given order p are cal- 
culated, one has no guarantee that all orbits of this order 
have been found. A numerical technique capable of cal- 
culating arbitrarily long periodic orbits to any desired ac- 
curacy was introduced in Ref. jlO| for the Henon map |Tl| l 
and later applied to a variety of other dynamical systems 
^2|-|l^|. Furthermore, this method provides a system- 
atic framework for the calculation of all periodic orbits 
of any given order p, in which each orbit is identified 
by a unique binary symbol sequence {s n }, n — I,..., p. 
The number of unstable periodic orbits of order p for the 
Henon map increases exponentially with p according to 
N(p) = 2 KoP , where Kq < 1 is the topological entropy. 



Therefore, the problem of finding all the periodic orbits 
of order p requires resources exponential in p. 

Series expansions over periodic orbits used for calcu- 
lations of dynamical averages are typically ordered ac- 
cording to the orbit length p [^]J|jl^,Q. Due to slow 
convergence and the fact that the number of orbits in- 
creases exponentially with p, these series often require 
a huge number of orbits |l5|,[l4]]. It was recently pro- 
posed (110 that using the cycle expansion framework 
for generic dynamical systems one can obtain better con- 
vergence by truncating the expansion according to the 
stability of the orbits [Q rather than their length p. This 
proposal is particularly useful since stability truncation 
does not require detailed understanding of the symbolic 
dynamics, it tends to preserve the shadowing properties 
and takes into account only the significant orbits of each 
length 0. 

In this Letter we show that a recently proposed tech- 
nique |nj provides a systematic framework for the cal- 
culation of the least unstable periodic orbits of any given 
order p. The resulting orbits are sorted according to 
their Lyapunov exponents starting with the least unsta- 
ble. The technique is highly flexible and can be applied in 
a straightforward manner to a great variety of discrete as 
well as continuous dynamical systems of any dimension. 
Therefore, it opens the way for employing the proposal 
of Refs. for a great variety of dynamical systems. 

We will first describe the method. Given a TV- 
dimensional chaotic dynamical system U : ri+\ = fifi) 
we generate a set of dynamical systems through the linear 
transformation: 



Sk ■ r i+ i = n + Ak[f(n) - n\ 



(1) 



where Afc are invertible TV x TV constant matrices which 
can be cast in the form = XCk with < A < 1. The 
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matrices Cfc are orthogonal with only one nonvanishing 
entry (±1) per row or column. The systems Sk are equiv- 
alent to the system U in the sense that there is a one to 
one correspondence between the fixed points of U and 
those of Sk ■ The important advantage of the representa- 
tion of Eq. (0) however is that, using a sufficiently small 
value for the parameter A, the fixed points of the trans- 
formed systems Sk become stable and therefore can easily 
be determined by an iterative process. Furthermore, the 
radius of convergence of this iterative algorithm turns out 
to be finite. The above procedure can be easily extented 
to the higher iterates f^ p '(r) of U [by replacing in Eq. 
(0) / with /( p )] allowing us to determine all the order 
p cycles of U . The parameter A is a key quantity here. 
For a given period p, it operates as a filter allowing the 
selective stabilization of only those unstable periodic or- 
bits, which possess Lyapunov exponents smaller than a 
critical value. Therefore, starting the search for unstable 
periodic orbits within a certain period p with a value of 
A = O(10 _1 ) and gradually lowering A we obtain the list 
of all unstable orbits of order p, starting with the least 
unstable orbit and sorted with increasing values of their 
Lyapunov exponents. 

Next, we focus on two dimensional systems. Denote 
the stability matrix of some periodic orbit of order p by 

(p) 

M, its matrix elements by rnf- where i,j — 1,2, and 
the stability eigenvalues of M by pf^ and p% . Without 



(p) 



> 



\p[ p) \. In Refs. 



loss of generality we assume \p{ 

|l9f a minimal set of matrices {Ck\ k = f,..,5 }, which 
is necessary and sufficient to achieve stabilization for any 
kind of hyperbolic fixed points was provided. The hyper- 
bolic fixed points with reflection (namely, fixed points for 
which at least one of the eigenvalues p{ , p 2 i s negative 



while \p (p) \ > I and \p {p> \ < 1) which satisfy pf> < 
become stable in the transformed system of Eq. p) if we 
use Ai = ACi with C\ = 1, where 1 is the unit 2x2 
matrix. After a little algebra a simple relation between 
the eigenvalues /4 2 °f the stabilized system Si and the 
eigenvalues p^ p \ m the original system U can be obtained: 



A«8 = 1 - A(l - pig). 



(2) 



The stability condition 1/4*2 1 < 1 leads to the critical 
value p c — 2/ A — 1 which represents an upper limit for 

the magnitude of the eigenvalues p{ 2 of the fixed points 
which are stabilized for a given A. This means that 

only those orbits for which \pf^\ < Pc become stable for 
the corresponding A. Therefore, varying A we can selec- 
tively extract those periodic orbits which possess stability 
eigenvalues less or equal to a given threshold value. 

The stabilization process for the case of hyperbolic 
fixed points with reflection and p\ > or without re- 
flection involves the matrices {Cfc | k = 2, .., 5} |L9). For 
these cases no exact monotonic relationship like Eq. (|^) 
can be derived. However, apart from exceptional cases 



(see below) , a selective stabilization procedure is possible 
similar to the case of the matrix C\\ the overall tendency 
is again that large values of A stabilize only the least un- 
stable periodic orbits and with decreasing value of A we 
get more and more of the increasingly unstable periodic 

orbits. To derive this let us consider without loss of gen- 
(p) (p) 

erality the case m{[ > m 22 . The matrix which has to 
be used for the stabilization process in this case is 



— I 
I 



(3) 



(p) 

The eigenvalues ^ of the transformed system S2 ex- 

pressed in terms of the eigenvalues p\ % of the original 
system are 



u [p) - 1 - - 



(4) 



(p) 
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,(p) 



Since typically p± ^> p 2 v> one 



where t — 

can use the approximation Tr(M) — (m^ + m^) = 
{p^ + p^) = p^ • Inserting this in the expressions for 



Pi P 2 and using [p^] 2 ^> Apf 1 p%' a detailed analysis of 

the possible cases leads to the statement p^ = amfj' 
where a is a factor of order unity. As can be seen 

(p) 

from p\ 2 the only exception to this situation occurs if 

m 22 ' l m< ii = ~ 1 + £ j where < e <C 1 which is cer- 
tainly rare in chaotic systems. We then have t = (3pf^ 
with j3 being a factor of 0(1). If the square root in 

Eq.(|) is real we have t 2 > 4(^ p) - 1)(1 - p ( 2 p) ) and 
it immediately follows that the relevant larger eigenvalue 
obeys /4 P = 1 — Xjpi /2 with 7 being a factor 
of O(l). If the square root is imaginary the real part 
of /i^ p \ which is responsible for the stabilization, obeys 
Re{p!^) = \ — \"fp^/2. Exceptional cases here are given 

by m 22 ^ / mfi = 1 — e where < e <C 1 . Although there is 
no strict monotonic ordering the above arguments clearly 
demonstrate that a monotonic ordering occurs to a good 
approximation. In addition we have performed a ran- 
dom matrix simulation for the original stability matri- 
ces respecting the above constraints due to hyperbolicity 
and calculated the distribution of the resulting prefac- 
tors 7 occurring in the eigenvalues of the transformed 
system. The results of this analysis confirm the above- 
obtained conclusions. The cases involving the other ma- 
trices {Cfc } can be treated analogously. 

To demonstrate the power of the above method we now 
apply it to the Henon map. To examine the spectrum of 
Lyapunov exponents we first calculated all the periodic 
orbits up to p = 23 for the Henon map ITT]] 
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-1 = a - x 



by n , 



yn-\-l — %n 



(•5) 



with the parameters a — 1.4 and b = 0.3, using the 
method described in [pTi|- The total number of orbits 
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obtained is 118,407 (including cyclic permutations and 
repetitions of lower cycles). For each one of these or- 
bits we calculated the Lyapunov exponent h = log(|p|)/p 
where p (in absolute value) is the largest eigenvalue of 
the matrix M = M p ■ . . . • M 2 • Mi, where 



M n = 



-2x n b 




1 



(6) 




are the Jacobian matrices of the map. The Lyapunov ex- 
ponents of all orbits of order p — 1, . . . , 23 as a function of 
p are shown in Fig. 1. We observe that in this attractor 
the two fixed points have the largest Lyapunov exponents 
and are thus the most unstable. The other Lyapunov ex- 
ponents form a band that becomes denser and broader 
as p increases. We also observe a small number of orbits 
with unusually small Lyapunov exponents. Such orbits 
appear for orders 13, 16, 18 and 20. 

To examine the dependence of Lyapunov exponents on 
the symbol sequence we plot the Lyapunov exponents for 
all the periodic orbits of order p = 20 vs. the sequential 
number of the orbit from to 2 20 — 1 (Fig. 2). The or- 
bits are ordered such that for every two adjacent orbits 
in the plot the symbol sequences are different in only one 
bit. This allows us to examine how the Lyapunov expo- 
nents change when one bit is switched in a long symbol 
sequence. To achieve this the symbol sequence \s n } for 
each orbit is considered as a Gray code sequence Q . A 
Gray to binary transformation {s„} — > {s' n } is then used 
and the decimal representation of {s' n } is given in the 
horizontal axis of Fig. 2. 
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FIG. 1. The Lyapunov exponents of all the periodic or- 
bits of order p = 1, . . . , 23 as a function of p for the Henon 
attractor at a — 1.4 and b = 0.3. 
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FIG. 2. The Lyapunov exponents for all the orbits of or- 
der p — 20 in the strange attractor at a = 1.4 and b = 0.3 
as a function of the (decimal representation of the) symbol 
sequence, converted into the Gray code. The vacant domains 
in the plot correspond to prunned periodic orbits. 

To stabilize the periodic orbits in the Henon attractor 
it turns out that one needs only two of the matrices 
namely C-y = 1 and C 3 = — C 2 . Using A- values in the 
range (0.05, 0.002) and a set of 500 starting points on 
the attractor we were able to find, for each of the periods 
p = 1, . . . , 23, within a few seconds of computation on 
a desktop workstation, the two least unstable periodic 
orbits. The results for p — 16, ... , 23 are presented in 
Table 1 which shows the period, the (cc, y) coordinates 
of one point of each orbit and its Lyapunov exponent 
hP = Qn\pM\)/p. 





Period 


x-coord. 


y-coord. 


Lowest Lyap. Exp. 


16 


1.414441 


0.525388 


0.261873 


16 


0.207904 


-1.293373 


0.259960 


17 


1.168372 


0.719909 


0.380900 


17 


0.956694 


0.775439 


0.379981 


18 


-1.255278 


1.588681 


0.285758 


18 


-1.256032 


1.588953 


0.284727 


19 


0.475407 


0.932648 


0.365918 


19 


0.608248 


0.683230 


0.365689 


20 


0.999632 


0.778785 


0.279102 


20 


0.687520 


-0.496536 


0.278732 


21 


1.433765 


-0.639919 


0.323827 


21 


1.018516 


0.377002 


0.323259 


22 


1.184577 


0.641833 


0.300455 


22 


0.641792 


0.655127 


0.300439 


23 


1.416001 


0.524053 


0.295087 


23 


1.596043 


-0.481690 


0.294950 



TABLE I. The two periodic orbits with the lowest Lya- 
punov exponents for orders p = 16, ... , 23, for the Henon 
map with a = 1.4 and b — 0.3. The x and y coordinates of 
one point of each orbit are shown as well as the Lyapunov 
exponent of the orbit. 
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FIG. 3. The Lyapunov exponent of each periodic orbits of 
period 16 for the Henon map (a = 1.4,6 = 0.3) is shown as 
a function of the critical value of the tuning parameter A be- 
low which the orbit is stabilized. Full circles represent orbits 
stabilized through Ci, while open circles represent orbits sta- 
bilized through C3. For Ci orbits, the Lyapunov exponents 
are strictly monotonic vs. A, while for C3 orbits only the 
general trend is monotonic with some deviations. 

To demonstrate how the periodic orbits are stabilized 
as the tuning parameter A is decreased we chose to 
present the results for the period p = 16. The total 
number of prime periodic orbits (namely, not including 
cyclic permutations and repetitions of smaller cycles) for 
p = 16 is 102. For a small value of A = 10~ 4 and 500 
starting points on the attractor we get all the 102 orbits 
stabilized. To examine the stabilization process we start 
with A = 0.05 and gradually lower it. This way, for each 
periodic orbit we identify the critical value of A below 
which the orbit is stabilized. In Fig. 3, we present for all 
orbits of period p = 16 the Lyapunov exponent of each 
orbit vs. the critical value of A below which this particu- 
lar orbit is stabilized. Orbits stabilized with the matrix 
C\ are shown in full circles, while orbits stabilized with 
C3 are shown in empty circles. We observe that for or- 
bits stabilized with C 1 the Lyapunov exponents increase 
monotonically as A is lowered. For orbits stabilized with 
C3 monotonicity is not strict, however the general trend 
is the same. The nearly monotonic tendency of the Lya- 
punov exponents vs. A demonstrates the suitability of 
our approach to determine, with varying A the set of 
least unstable periodic orbits within a given period. 

We observe that the spectrum of Lyapunov exponents 
(Fig. 2) exhibits a rough landscape, which resembles the 
energy spectrum which typically appears in hard mini- 
mization problems. Well known problems of this type are 
finding spin glass ground states |]2l| and protein folding. 
The use of combinatorial techniques for these problems is 
infeasible since the computational resources required are 
exponential in the input size. However, probabilistic al- 
gorithms such as simulated annealing |22] can provide ap- 



proximate results, which for many applications are prac- 
tically sufficient. Our results may provide new insight 
about hard minimization problems in general. The issue 
is whether one can artificially destabilize all the meta- 
stable states above some externally tuned energy E and 
this way accelerate the convergence into the low-lying 
states. 

In summary, we have shown that using the method 
proposed in Refs. jl9| one can obtain the periodic or- 
bits of any given order p sorted according to their Lya- 
punov exponents starting with the least unstable one. 
The method can be applied to a great variety of discrete 
as well as continuous dynamical systems of any dimen- 
sion. Having the periodic orbits sorted in increasing order 
of their Lyapunov exponents is highly useful in light of 
the recent proposal [ ^6|Jl7| that in cycle expansion calcu- 
lations for generic dynamical systems better convergence 
can be obtained by truncating the expansion according 
to the stability of the orbits rather than their length p. In 
particular, stability truncation does not require detailed 
understanding of the symbolic dynamics, it tends to pre- 
serve the shadowing properties and takes into account 
only the significant orbits of each length, leaving out an 
exponential number of insignificant orbits. 
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